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Recently, ultrastable glasses have been created through vapor deposition. Subsequently, computer 
simulation algorithms have been proposed that mimic the vapor deposition process and result in 
simulated glasses with increased stability. In addition, random pinning has been used to generate 
very stable glassy configurations without the need for lengthy annealing or special algorithms in¬ 
spired by vapor deposition. Kinetic and mechanical stability of experimental ultrastable glasses is 
compared to those of experimental glasses formed by cooling. We provide the basis for a similar 
comparison for simulated stable glasses: we analyze the kinetic and mechanical stability of simu¬ 
lated glasses formed by cooling at a constant rate by examining the transformation time to a liquid 
upon rapid re-heating, the inherent structure energies, and the shear modulus. The kinetic and 
structural stability increases slowly with decreasing cooling rate. The methods outlined here can be 
used to assess kinetic and mechanical stability of simulated glasses generated by using specialized 
algorithms. 

PACS numbers: 64.70.Q-,81.05.Kf 


Vapor deposition on a substrate held at around 85% of 
the glass transition temperature is used to create glasses 
that have a higher kinetic stability [Tj and larger elastic 
moduli [2] than glasses created by cooling from a liq¬ 
uid. Due to these and other desirable properties of these 
so-called ultrastable glasses, it is of great interest to un¬ 
derstand the differences between glasses created through 
vapor deposition and by cooling from a liquid. Since va¬ 
por deposited glasses are created one layer at a time, it 
is suspected that enhanced mobility at the surface allows 
the molecules to find more stable configurations. 

Discovery of ultrastable glasses in the lab inspired sev¬ 
eral computer simulation studies. The goal of these stud¬ 
ies has been two-fold. First, modeling the vapor depo¬ 
sition process can provide insight into the origin of the 
stability of vapor deposited glasses. Second, the experi¬ 
mental procedure can be used as an inspiration for the de¬ 
velopment of computer simulation algorithms that could 
generate very stable simulated glasses without the need 
for lengthy annealing. Leonard and Harrowell [3] mod¬ 
eled ultrastable glasses using a three-spin facilitated Ising 
model. Their study provided evidence that enhanced sur¬ 
face mobility during vapor deposition results in improved 
kinetic stability, but the simplicity of their model pre¬ 
vented them from studying other properties of ultrastable 
glasses, such as the shear modulus. On the other end of 
the complexity spectrum, Singh et al. [3] modeled vapor 
deposition of the organic molecule trehalose by introduc¬ 
ing 1-5 “hot” molecules above a substrate and then (arti¬ 
ficially) cooling them to the substrate temperature. This 
procedure resulted in a simulated glass with a higher ki¬ 
netic stability, but also pronounced structural anisotropy 
in the direction perpendicular to the substrate. 

More recent studies used well-known glass-forming bi¬ 
nary Lennard-Jones mixtures [5]. First, Singh, Ediger, 


and de Pablo [6] and Lyubimov, Ediger, and de Pablo 
[7] modeled the vapor deposition of a binary mixture of 
Lennard-Jones particles to simulate the creation of an ul¬ 
trastable glass. They found an enhanced kinetic stability 
compared to the glass obtained by cooling from a liquid 
at a constant rate, but the enhanced stability was less 
pronounced than in experiments. Furthermore, the simu¬ 
lated vapor deposition resulted in glasses with structural 
anisotropy, a higher density, and non-uniform composi¬ 
tion compared to the simulated glasses created by cool¬ 
ing from a liquid. Very recently, Hocky et al. |5] created 
two dimensional binary Lennard-Jones mixture glasses 
by randomly pinning particles. Glasses formed by ran¬ 
dom pinning are, by construction, isotropic and uniform 
(on average). Hocky et al. found that these glasses have 
increased kinetic stability but the presence of pinned par¬ 
ticles prevented them from investigating whether these 
glasses also have a larger shear modulus. 

The above mentioned simulational studies focused on 
the properties of stable glasses created by using different 
specialized algorithms. However, it is difficult to assess 
the stability of these glasses without knowledge of the 
stability of simulated glasses generated in a conventional 
way, i.e. by cooling at a constant rate. We emphasize 
that simulated glasses generated by cooling are isotropic 
and uniform. Thus, studying these glasses may also shed 
some light on the importance of the anisotropy and com¬ 
positional non-uniformity for the glass stability. Here we 
create glasses by cooling from the supercooled liquid and 
we assess their kinetic and mechanical stability by exam¬ 
ining their kinetic stability against reheating, properties 
of their potential energy landscape, and their shear mod¬ 
ulus. 

We simulated the Kob-Andersen (KA) [5] 80:20 
binary Lennard-Jones mixture in three dimensions. 
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The particles interact via the potential V a p (r) = 

4 e a p (-^-) 12 — j and the parameters are cbb = 


0.56,4,4, CAB = 1-56AA, cr BB = 0.88ctaa, and (TAB = 
0.8ctaa- The majority species are of type A and both 
masses are equal. We cut the potential at 2.5 a a p. The 
results are presented in reduced units with cfaa, caa/^b, 
and \Jm,Acr 2 AA /caa being the units for length, temper¬ 
ature, and time, respectively. We simulated N = 8000 
particles at a number density N/V = 1.20 using periodic 
boundary conditions. 


We ran NVT simulations using a Nose-Hoover thermo¬ 
stat. We started with an equilibrated supercooled liquid 
at the temperature T = 0.5 (for this system the onset 
temperature for the slow dynamics is T a ss 1.0 and the 
mode-coupling temperature is T c = 0.435). We ran 80 
cooling runs at cooling rates T = AT/At of 3.33 x 10~ n 
where n = 3, 4, 5, 6, and 7. Each run started at an 
equilibrium configuration at T = 0.5 and was cooled to 
T = 0.3. We also ran 4 cooling runs at a cooling rate 
of 3.33 x 10~ 8 . For T > 3.33 x 10 -8 , we started from 
the final configuration of each of the 80 cooling runs and 
annealed the system for at least 100 times r Q determined 
from T = 0.5 (we define T a through the average over¬ 
lap function discussed later). For the cooling rate of 
3.33 x 10” 8 , we performed 5 different annealing runs for 
each of the 4 configurations, each beginning with a dif¬ 
ferent initial random velocity. In addition, we heated 
the final configurations of the cooling runs by ramping 
up the temperature from T = 0.3 to T = 0.5 at a con¬ 
stant rate over a time t = 10 (we tried increasing the 
temperature instantaneously from T = 0.3 to T = 0.5 
but found that this resulted in large, non-physical, os¬ 
cillations of the kinetic and potential energies). After 
ramping up the temperature to T = 0.5, we continued 
constant temperature simulations (at T = 0.5) until the 
mean-square displacement started growing linearly with 
time. We refer to this ramping up of temperature and 
the subsequent run at T = 0.5 as a heating trajectory. 
For T > 3.33 x 10~ 8 , we ran 80 heating trajectories, in 
each case starting from the final configuration of a differ¬ 
ent cooling run. For T = 3.33 x 10~ 8 , we ran the heating 
trajectories 15 times using different initial random veloc¬ 
ities for each of the configurations obtained from the 4 
cooling runs. Finally, we also performed an equilibrium 
run at T = 0.5, for at least 100r Q . We used HOOMD- 
blue um for equilibration at T = 0.5 and the cooling 
simulations. The heating trajectories, the T = 0.3 runs, 
and the single run at T = 0.5 used LAMMPS [HUH] run 
on a GPU jlBKIB] . 

We start with an assessment of the kinetic stability of 
the simulated glasses. First, we we examine the time it 
takes for the particles to rearrange after the sudden in¬ 
crease in temperature, i.e. along the heating trajectories. 
We quantify these rearrangements through the average 
overlap function q s (t,tw), which measures the probabil- 



FIG. 1: Average overlap function q s (t,tw = 0) monitored 
during the heating trajectories (solid lines) compared to the 
average overlap function for the equilibrium liquid at T = 0.5 
(dashed line). Solid lines correspond to glasses obtained with 
cooling rates T = 3.33 x 10 -71 , where n = 3, 4, 5, 6, 7, 
and 8 from left to right. The inset shows the mean square 
displacement versus time for the same cooling rates (solid 
lines) and for the equilibrium liquid at T = 0.5 (dashed line). 


ity that a particle moved over distance a during the time 
between t w and t + t w , where t w is the waiting time. 



where qmi/jtw') — 0 (u |r m (t T t^,) t*m(^i(j)|), 0 is 

Heaviside’s step function, and r m (t) is the position of 
a particle to at a time t. t w , the waiting time, is mea¬ 
sured from the start of the heating trajectory. We chose 
a = 0.25 to be consistent with previous work [TB]. We 
define r a as when q s (r a ,tw) = e _1 for an equilibrium 
system (note that for an equilibrium system q s (t,tw) 
does not depend on the waiting time, i.e. it is time- 
translationally invariant). 

Fig. [I] shows average overlap function q s (t,tw = 0) ob¬ 
tained while heating glasses prepared at different cooling 
rates (solid lines) and obtained from the equilibrium run 
at T = 0.5 (dashed line). The feature at t = 10 is a 
consequence of changing the thermostat from a ramp-up 
of the temperature at a constant rate from T = 0.3 to 
T = 0.5 to maintaining constant temperature T = 0.5. 
For t > 10, q s (t,t w = 0) exhibits prolonged plateaus for 
the three slowest cooling rates followed by a very rapid 
decay at later times (the glass generated using the slowest 
cooling rate,3.33 x 10 —8 , relaxes faster than for ballistic 
motion). With decreasing cooling rate both the plateau 
height and its extent increase. This indicates that with 
decreasing cooling rate the cage diameter decreases and 
it takes longer for the particles to move out of their cages. 
This is the first indication of increasing kinetic stability 
against re-heating. 

The inset in Fig. [l] shows the mean square displace- 
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FIG. 2: Dependence of stability ratio S = ttrans/T a on cooling 
rate T. The inset shows the average overlap function versus 
time for the cooling rate T = 3.33 x 10 -6 and waiting times, 
from top to bottom, 0, 10, 500, 1000, 3000, 4250, 5000. 

ment ( 6r 2 (t,t w )) = N _1 (Y, n l r n{t + t w ) - r n (t w )] 2 ^ 
for a waiting time t w = 0 obtained while heating glasses 
prepared at different cooling rates (solid lines) and ob¬ 
tained from the equilibrium run at T = 0.5 (dashed line). 
We find that ( Sr 2 (t,t w )) encodes similar information as 
Qs{t,t w =0). In particular, we note that with decreas¬ 
ing cooling rate the time required for the mean-square 
displacement to reach the equilibrium curve increases. 
This is another indication of increasing kinetic stability 
against re-heating. 

To quantify kinetic stability we follow earlier experi¬ 
mental m and simulational [8] studies and we evaluate 
the transformation time t tra ns and the stability ratio S. 
The transformation time is defined as the time it takes 
for the liquid obtained by heating the glass to return 
to equilibrium. In the experimental study of Ref. m 
the transformation time was obtained by monitoring the 
time evolution of the dielectric response. We follow the 
simulational study of Ref. [8] and obtain the transforma¬ 
tion time from monitoring the waiting time dependence 
of the average overlap function. Specifically, we calcu¬ 
late q s (t, t w ) for a range of t w (see the inset to Fig. [2] for 
q s {t,t w ) for various waiting times for T = 3.33 x 10” 6 ) 
and we define the waiting time-dependent relaxation time 
r s (t w ) through the equation q s (T s ,t w ) = e -1 . We define 
the transformation time, ttrans , as the waiting time such 
that T s (t u ,) = T a , where r a is the equilibrium relaxation 
time. Finally, we follow Refs. 0113 and we define a 
stability ratio S = t tra ns/ T a- 

In Fig. [2] we show the cooling rate dependence of the 
stability ratio. At larger cooling rates S increases rather 
quickly with decreasing T. However, for cooling rates less 
than 3.33 x 10~ 5 we get a slower increase of the stability 
ratio with cooling rate. For the slowest cooling rates, an 
order of magnitude decrease in the cooling rate results 
in a doubling of the stability ratio. Thus, to achieve sta¬ 



FIG. 3: Average potential per particle, (U), as a function of 
cooling rate, T. The inset shows average inherent structure 
energy per particle, Eis, as a function of cooling rate. The 
calculations were done at cooling rates 3.33 x 10“ n , with n = 
8, 7, 6, 5,4, 3, the lines are to guide the eye. 


bility ratios of the most stable simulated glasses, which 
approach approximately 400 [8], we would need to de¬ 
crease the cooling rate by about 3 decades (to achieve 
stability ratios of experimental ultrastable glasses, which 
approach approximately 10 3 5 [17] we would need to de¬ 
crease the cooling rate by about 22 decades). We con¬ 
clude that even simulated stable glasses would be impos¬ 
sible to obtain by cooling, at least with the computational 
resources available to us at present. 

Next, we turn to an assessment of the mechanical sta¬ 
bility of the simulated glasses. The simplest, albeit in¬ 
direct, measure of the mechanical stability is provided 
by the inherent structure energy [18] , since in the en¬ 
ergy landscape picture of the glass transition, more sta¬ 
ble glasses are in a lower energy basin of the landscape 
and have to overcome a larger energy barrier to deform. 
We show the dependence of the average inherent struc¬ 
ture energy (obtained using the Fire algorithm [19 in 
HOOMD-blue) at T = 0.3 on the cooling rate used to 
reach this temperature in the inset in Fig. [3] We see that 
the average inherent structure energy decreases with de¬ 
creasing cooling rate. In Fig. [3] we show the dependence 
of the average potential energy per particle at T = 0.3 on 
the cooling rate. It also decreases with decreasing cool¬ 
ing rate. We note that for their vapor deposited glass at 
T = 0.3, Ediger and de Pablo got an inherent structure 
energy of -8.35 [7] and a potential energy of -7.8 (esti¬ 
mated from [ 7 ]). These values are close to our values of 
-8.30 and -7.86. We note, however, that the vapor de¬ 
position inspired algorithm used in Ref. [7 resulted in 
glasses that were non-uniform in contrast to our simu¬ 
lated glasses, and thus, a literal comparison of the two 
studies is impossible. 

A more direct measure of mechanical stability is pro¬ 
vided by elastic constants. Ultrastable glasses prepared 
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in experiments [2] and stable glasses obtained in simula¬ 
tions using vapor deposition-inspired algorithms [6j have 
larger elastic constants than glasses prepared by ordinary 
cooling. We examine here the cooling rate dependence of 
the shear modulus of our simulated glasses. To this end 
we use an approach recently introduced by two of us [20]. 
This method is based on the relationship between corre¬ 
lations of transverse particle displacements in the glass 
and the glass’s shear modulus. The final result is that 
the shear modulus of a glass can be measured by study¬ 
ing the small wave-vector q behavior at long times of the 
four-point structure factor 

SM;t) = 1 /j2g[Sr n (t)]g* [6r m (t)] e ^'1^(°)-^(°)1 \ , 

\n,m / 

( 2 ) 

where the weighting function g[5r n (t)] = 5r^{t) = 
i r n(t) ~ r n (0))' L and Sr^(t) is the component of 
the displacement of particle n perpendicular to the 
wave-vector q, Sr^(t) ■ q = 0. As argued in 

Ref. [20], the shear modulus for a glass can be ob¬ 
tained from the small q , long t behavior of S 4 (q;t), 
/i = lim^o linii-Kxj 2fcgTp [q 2 S±(q] t)] . In practice 

one can examine the small wave-vector behavior of 
2kgTp [q 2 S^(q; t)] for times longer than the initial 
transient dynamics but shorter than time scale for the 
aging process of a glass. Also, in order to get good 
statistics one either needs many more independent con¬ 
figurations at the end of each cooling process than the 80 
that we generated or one needs to continue running at 
T = 0.3 and average over different time origins. In order 
to measure properties of the glass as prepared using a 
specific cooling rate, one can only average over different 
time origins before the start of the aging process. For 
faster cooling rates the aging starts very quickly and the 
range of available time origins is rather small. Fortu¬ 
nately, the present method to calculate the shear mod¬ 
ulus does not require very long trajectories [2D] - We 
show 2 kBTp[q 2 S 4 : (q- 1 t)\ versus q at a time t = 100 
for T = 0.3 in the inset to Fig. [4] We averaged over time 
origins for t < 100 for the 3 faster cooling rates, and for 
t < 400 for the slowest cooling rate. 

In Fig. [4] we show the cooling rate dependence of the 
shear modulus. We find that the modulus increases with 
decreasing cooling rate, although the rate of increase 
seems to be decreasing with decreasing T. At the slowest 
cooling rates an order of magnitude decrease of the cool¬ 
ing rate results in an increase of the modulus by 4.6%. 
We note that Ashwin, Bouchbinder, and Procaccia [21] 
also investigated the cooling rate dependence of the shear 
modulus. They used a different system, in two spatial di¬ 
mensions and thus their results cannot be literally com¬ 
pared to ours. However, we do note that in a similar 
range of cooling rates our shear modulus increases by 
15% and theirs increases by 12%. Interestingly, these 



FIG. 4: Dependence of the shear modulus, p, on the cool¬ 
ing rate T for glasses at T = 0.3. The inset shows 
2 k B Tp [q 2 S4q-,t)] versus q plotted at a time t = 100 for 
T = 0.3, for the four slowest cooling rates: 3.33 x 10 -n where 
n = 5,6, 7,8, from bottom to top. 


relative increases of the shear modulus achieved in our 
and Ashwin et al. simulational studies are comparable 
to the relative differences between the shear modulus of 
experimental ultrastable glasses and experiemntal glasses 
generated by cooling (19% for indomethacin and 15% for 
tris-naphtylbenzene ED- It is unclear to us whether this 
fact has some deeper meaning. 

In summary, we showed here that by decreasing the 
cooling rate in simulations, we generated simulated 
glasses that have increasing kinetic and mechanical sta¬ 
bility. The highest kinetic stability that we achieved is 
about 66. Decreasing the cooling rate by one decade, 
which would result in the kinetic stability of about 130, 
would require approximately 1 year on a single GPU us¬ 
ing HOOMD-blue per one cooling trajectory. This re¬ 
sult establishes a baseline against which specialized al¬ 
gorithms that are inspired by vapor deposition or use 
random pinning should be measured. 

We believe that the kinetic stability, quantified by the 
stability ratio S = t tra ns/Ta is the best quantity to com¬ 
pare stability of glasses generated using different proto¬ 
cols and algorithms. We note that the average potential 
energy is very sensitive to the average density, compo¬ 
sitional nonuniformities and possibly to anisotropy. In 
fact, our glasses generated by cooling at constant rate 
have lower average potential energy than stable glasses 
generated by Singh et al. [6j using a vapor deposition 
inspired algorithm. We note the average inherent struc¬ 
ture energy seems to correlate better with the stability 
but most likely it is also sensitive to the density and 
compositional nonuniformities. Finally, we note that the 
relative change of the shear modulus of simulated glasses 
formed by cooling at different cooling rates is compara¬ 
ble to the relative difference between the shear modulus 
of experimental glasses created through vapor deposition 
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and by cooling. This surprising finding, which deserves 
further investigation, suggests that that the kinetic and 
mechanical stability are not necessarily correlated. 

The methods presented here can be used to assess ki¬ 
netic and mechanical stability of simulated glasses gen¬ 
erated by using specialized algorithms. 

We gratefully acknowledge the support of NSF grant 
CHE 1213401. 
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